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ABSTRACT 

The multi-planetary system HD 128311 hosts at least two planets. Its dynamical formation 
history has been studied extensively in the literature. We reanalyse the latest radial velocity 
data for this system with the affine-invariant Markov chain Monte Carlo sampler EMCEE. Us¬ 
ing the high order integrator IAS 15, we perform a fully dynamical fit, allowing the planets to 
interact during the sampling process. A stability analysis using the MEGNO indicator reveals 
that the system is located in a stable island of the parameter space. In contrast to a previous 
study, we find that the system is locked in a 2:1 mean motion resonance. The resonant angle 
(fi is librating with a libration amplitude of approximately 37°. The existence of mean mo¬ 
tion resonances has important implication for planet formation theories. Our results confirm 
predictions of models involving planet migration and stochastic forces. 

Key words: methods: numerical — gravitation — planets and satellites: dynamical evolution 
and stability — stars: individual (HD 128311) 


1 THE PLANETARY SYSTEM HD128311 


The first planet around HD 128311 was discovered by [Butler et al. 

S . A second planet was found two years later ( jVogt et al. 

Both planets are most likely gas giants with minimum 
masses of 1.8 Mjup and 3.2 Mjup. A third signal has been discov¬ 
ered by jMcArthur et al.| ( |2014| ) but its planetary nature has yet to be 
confirmed. Soon after its discovery HD 128311 began to emerge 
as a test bed of planet formation scenarios within the scientific 
community. A large number of groups studied the formation and 
evolution of this system with particular focus on the system’s dy- 


namical history ([Beauge et al.|2006|[Quillen|2006|[ 

Sandor & Kley 

2006 

[Sandor et al. [20071 [Michtchenko et al.|20081[Meschiari et al. 

2009 

[Crida et al.[[2008| [Lee & Thommes[[2009| 

Voyatzis et al. 

2009 

[Barnes & Greenberg[[2006| [Raymond et al. 

[20081 [Rein & 

Papa 

oizou[2009 [Lecoanet et al.[2009[[Gavon-Markt & Bois [20091 

Gozdziewski & Konacki[[2006 [Erdi et al.[[2007 Giuppone et al. 

2010[>. 


The most important dynamical feature of HD 128311 is its 
proximity to a 2:1 mean motion resonance (MMR). Early observa¬ 
tions supported the idea that HD 128311 is in resonance, although 
some of the orbital solutions were dynamically unstable on short 
timescales ( Vogt et al.|2005| ). The most recent work on HD128311 
by [McArthur et al. ( [2014| ) includes new radial velocity data, a re¬ 
calibration of older data sets and astrometric and photometric con¬ 
straints. The authors found that the planets in their best fit model 
are not in a MMR. Whether this system is in a MMR or not is an 
important constraint for planet formation scenarios. For example 


a system in a MMR favours the idea that giant planets migrated 
when they were still embedded in a protoplanetary disc ( [Lee &[ 
[Peale|20()^[Rein & Papaloizou|2009[ ). 

In this letter we reanalyse the combined radial velocity (RV) 
datasets including the recalibrated data from the Hobby-Eberly 
Telescope (HET) and the Lick Observatory. Our RV data is thereby 
equivalent to that of [McArthur et al.[ ( |2014| ) but we do not take 
into account the astrometric observations in our model which only 
constrain the system’s inclination. In contrast to previous stud¬ 
ies, we use a fully dynamical model, allowing planets to interact, 
rather than being on fixed Keplerian orbits. We use a use a mod¬ 
ern Markov chain Monte Carlo sampler to explore the high dimen¬ 
sional parameters space and converge on a new set of orbital pa¬ 
rameters. The details are described in Section We then perform 
long term orbit integrations to study the stability in Sec tion [^ us- 
ing the fast chaos indictor MEGNO ( [Cincotta & Sim6[[2000| ). In 
Sectionj^we conclude by discussing the resonant nature of the sys¬ 
tem and the implications for planet formation scenarios involving 
stochastic migration. 


2 MARKOV CHAIN MONTE CARLO 

We use a Markov chain Monte Carlo method to fit the observed 
radial velocity datasets. In comparison to the fit performed by 
[McArthur et al.[ ( [20T4] ) we use a fully dynamical two planet model. 
For this purpose, we couple the high order IAS 15 ( [Rein & Spiegel[ 
[2015j ) integrator which is available within the REBOUND package 
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Figure 1. Radial velocity data and model. The time is relative to JD2450983.83. Top panel: RV data points, MCMC mean (red), MCMC samples (gray). 
Bottom panel: residual errors not including jitter 5 . 


( [Rein & Liu|2012| | to the EMCEE package. EMCEE is an open-source, 
parallelized affine-invariant Marcov chain Monte Carlo sampler 
written in python ( [Foreman-Mackey et al.|2013] ). 

We assume co-planarity of the planets but allow the orbital 
plane to be inclined from the line of sight by an angle i. The cen¬ 
tral object has a fixed mass of - 0.828Mo. We have four or¬ 
bital elements and one mass parameter, m sin(/), for each planet. 
The orbital parameters are the period P, the eccentricity e, the ar¬ 
gument of periapsis co and the mean anomaly M. For the sampling 
process, we perform a coordinate transformation io h - e sin(du) 
and k = e cos(du) variables. This allows us to avoid the coordinate 
singularity at ^ = 0 and speed up convergence (for a discussion see 
|Hou et aL|2012| l. We use a Jacobi coordinates, i.e. the outer planet’s 
orbital parameters are given with respect to the centre of mass of 
the star and the inner planet. Building up on the work of |McArthur| 
|et al.| ( [MT^ , we further include an individual offset y and jitter pa¬ 
rameter V for each RV dataset (four for Lick, one for HET). This 
allows the model to account for variations in the offsets for differ¬ 
ent instruments as well as variations in the Lick detectors between 
observing runs. The total number of free parameters is thus 21. 

A total of 500 walkers are evolved for several thousand gener¬ 
ations. Both angles co and M as well as the eccentricity e for each 
planet have a uniform priors. We follow [Hou et al.| ( [2012| ) and use 
uninformative Jefferys priors for the period, the planet mass as well 
as the jitter parameters. We initialize the walkers with mass and pe¬ 
riod parameters that are roughly those of previous results to speed 
up convergence. Our experiments showed that the precise details on 
how the walkers are initialized do not matter. Finally, we perform 
a simple declustering algorithm after 1000 generations by remov¬ 
ing those samples from the ensemble that have a likelihood signifi¬ 
cantly smaller than the best samples ( |Hou et al.|2012| ). 

After the MCMC has converged and has been declustered, we 
sample it over 1250 generations and calculate the mean of all pa¬ 
rameters as well as the two sigma confidence intervals. The results 
are listed in Table In Figure we plot the observed radial ve¬ 
locity datapoints with the offset adjusted according to our model. 
The red line corresponds to the mean solution of TableThe gray 


lines correspond to 50 randomly drawn samples from the MCMC 
posterior distribution. 


3 LONG TERM STABILITY 


We focus on parameters close to those of the converged MCMC 
samples and run a total of 57 600 realisations of HD 128311 for 
^max = 10^ years to map out the structure of the phase space. A 
smaller subset was integrated for ^max = 10^ years, but no qual¬ 
itative differences could be observed. As in the previous section 
we use REBOUND and the high order IAS 15 integrator for these long 
term integrations. We declare a system unstable if at least one of 
the planets gets ejected, if the planets have a close encounter, or 
if the semi-major axis of at least one planet changes by more than 
50%. We refer to this definition as Lagrange stability and call the 
time until instability Lagrange timescale. 

For systems that are Lagrange stable, we measure the Mean 
Exponential Growth of Nearby Orbits, MEGNO (|Cincotta & Sim6| 
12000] ). MEGNO is a fast chaos indicator similar to the tradi¬ 
tional Lyapunov exponent, but gives more useful results on shorter 
timescales. We compute the MEGNO value {Y) by integrating the 
variational equations { [Mikkola & Innanen|1999j ) using IASI5. Eor a 
de finition of (7) and a detail ed d iscussion of the MEGNO ind icator 
see Cincotta & Simo < |2000| ) and Gozdziewski et al. ( 200l|^ 


In Figure we show a slice of the parameter space in the 
plane of the outer planet’s period and eccentricity, P 2 and 62 . All 
other orbital parameters are those listed in Table All simula¬ 
tions in red regions are Lagrange unstable. The shade of red in¬ 
dicates the Lagrange timescale. Simulations in blue regions remain 
stable for the entire integration. The shade of blue corresponds to 
the MEGNO (F). A value of (F) = 2 (dark blue) indicates a stable 
quasi-periodic orbit, whereas a value larger than 2 indicates chaotic 
behaviour. One can see in the figure that seemingly stable systems 


^ Note that there is a typo in the denominator of S in Equation 5 of 
[Gozdziewski et al.H2001| . 
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Figure 2. Lagrange stability plot in the vicinity of the MCMC solution. The axes correspond to the period and eccentricity of the outer planet. Red regions are 
Lagrange unstable. Blue regions are stable for the duration of the integration. Dark blue regions are stable quasi-periodic solutions according to the MEGNO 
indicator. The white cross shows the nominal mean orbital parameters listed in Table[T] The black dots show samples of the MCMC posterior. 


which are close to the stability boundary are in fact chaotic and thus 
might become unstable over longer timescales. 

The white cross in Figurej^shows the mean orbital parameters 
of the MCMC sample. The black dots show the individual MCMC 
samples from the posterior distribution. Both the mean solution and 
all MCM samples are located well within a stable island. Note that 
the eccentricity 62 is not well constrained. However, the observed 
RV data is inconsistent with the assumption of a circular orbit. 

The observational data is currently not good enough to con¬ 
strain the mutual inclination of the two planets. To test the effects 
of mutually inclined orbits, we ran an additional 14 400 simulations 
with the same initial conditions as before but perturbed each system 
with a random mutual inclination of ~ 2°. Our results indicate that 
the system’s stability is not significantly affected by a small amount 
of mutual inclination over a million year timescale. 


4 MEAN MOTION RESONANCE AND CONCLUSIONS 

The period ratio in our MCMC sample, P 2 IP 1 , is within 2% of a 2:1 
mean motion resonance. The period ratio alone is not a sufficient 
criteria for a MMR. To test whether the system is in a MMR or 
not we therefore monitor the two resonant angles (fi - Ai - 2 A 2 + 
CDi and (P 2 = Ai - 2 A 2 + (D 2 , where A is the mean longitude. In 
all of our MCMC samples, ipi is librating around 0°. The libration 
amplitudes range from close to 0° to up to 90° with a mean of 37°. 
The resonant angle ip 2 is librating in the majority but not all of the 
samples. The mean libration amplitude for ip 2 is 76.5°. Our findings 
differ qualitatively from those of |Mc Arthur et al.| ( |2014| ) who found 
that the planets are not in a MMR. 


The existence of a MMR is an important observable in many 
planet formation scenarios. A system in a MMR supports the idea 
that giant planets migrated into their current position while they 
were still embedded in a protoplanetary disc ( |Lee & Peale|2002| . 
Furthermore, [Rein & Papaloizou] ( |2009| ) predicted that if the sys¬ 
tem did undergo a phase of stochastic migration, then ipi should 
librate at moderate amplitudes whereas ^2 should be close to the 
separatrix of libration. Our results are in perfect agreement with 
this prediction, thus supporting the idea that both planet migration 
and stochastic forces occurred during the system’s evolution. 
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Table 1. Model parameters from MCMC sample. 

Parameter 


Value and 2-cr conhdence interval 

epoch 

JD 

2450983.83 (fixed) 

stellar mass 

M, 

0.828 Mo (fixed) 

inclination 

i 

63.8!|;^° 

Planet 1 

minimum mass 

m\ sin(0 

1.83!«:15 Mjup 

period 

Pi 

460. g days 

eccentricity 

ei 

0 30^** **^ 

-0.04 

argument of periapsis 

OJl 

-76.2+®;^ ” 

mean anomaly 

Ml 

OCQ 9 + 11.9 o 
z,jy.z,i26 

Planet 2 

minimum mass 

m 2 sin(0 

8 M- 

period 

Pi 

910.7^g Q days 

eccentricity 

ei 

0 12+008 
-0.06 

argument of periapsis 

0^2 

-19 7+23-2 o 

mean anomaly 

M2 

184.2+20-0 “ 

Lick Radial Velocities starting JD 2450983 

offseP 

n 

1.2+l>m/s 

jitter*^ 

>^1 

1.5^j 4 m/s 

Lick Radial Velocities starting JD 2451409 

offseP 

72 

7.9!ll>/s 

jitter*' 

>^2 

5.2+10/ m/s 

Lick Radial Velocities starting JD 2452333 

offseP 

73 

11.2+2‘;0m/s 

jitter*' 

>^3 

6.4+^y m/s 

Lick Radial Velocities starting JD 2452515 

offseP 

74 


jitter*' 

54 

2 . 8 +fo m/s 

HET Relative Radial Velocities 


offset^ 75 0 - 6 ^} 4 m/s 

jitter*^ >^5 ^-^-2 0 


Notes; 

^ Offsets are relative to those found by [McArthur et al.| ( |2014) and are all 
consistent with zero. 

^ Jitter is added on top of the reported RV errors which do not include stellar 
variability and additional instrumental noise. 
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